Comparative classical and "ab initio" Molecular 
Dynamics study of molten and glassy germanium 
dioxide 

M Hawlitzky 1 , J Horbach 1 2 , S Ispas 3 , M Krack 4 and K 
Binder 1 

-^nstitut fiir Physik, Johannes Gutenberg-Universitat Mainz, Staudinger Weg 7, 
55099 Mainz, Germany 

2 Institut fiir Materialphysik im Weltraum, Deutsches Zentrum fiir Luft- und 
Raumfahrt (DLR), 51170 K61n, Germany 

3 Laboratoire des Colloi'des, Verres et Nanomateriaux, Universite Montpellier II 
and CNRS UMR 5587, 34095 Montpellier, France 

4 Computational Science, Dept. of Chemistry and Applied Biosciences, ETH 
Zurich, USI Campus, Via Giuseppe Buffi 13, 6900 Lugano, Switzerland 

E-mail: kurt.binder@uni-mainz.de and juergen.horbach@dlr.de 

Abstract. A Molecular Dynamics (MD) study of static and dynamic properties 
of molten and glassy germanium dioxide is presented. The interactions between 
the atoms are modelled by the classical pair potential proposed by Oeffner and 
Elliott (OE) [Oeffner R D and Elliott S R 1998 Phys. Rev. B 58 14791]. We 
compare our results to experiments and previous simulations. In addition, an 
"ab initio" method, the so-called Car-Parrinello Molecular Dynamics (CPMD), 
is applied to check the accuracy of the structural properties, as obtained by the 
classical MD simulations with the OE potential. As in a similar study for Si02, 
the structure predicted by CPMD is only slightly softer than that resulting from 
the classical MD. In contrast to earlier simulations, both the static structure 
and dynamic properties are in very good agreement with pertinent experimental 
data. MD simulations with the OE potential are also used to study the relaxation 
dynamics. As previously found for Si02, for high temperatures the dynamics of 
molten Ge02 is compatible with a description in terms of mode coupling theory. 



PACS numbers: PACS numbers: 61.20.Lc, 61.20.Ja, 02.70.Ns, 64.70.Pf 

1. Introduction 

Understanding the structure and dynamics of glassforming fluids and the nature of 
the glass transition is one of the most challenging unsolved problems of the physics of 
condensed matter Q] [2l [3j [H El [7] . One of the most debated issues is the question 
to which extent the glass transition is a universal phenomenon; i.e. it is debated 
whether the mechanisms causing the dramatic slowing down in undercooled fluids 
when the glass transition is approached are basically the same in all glassforming 
materials, or whether qualitatively different classes of glass transitions exist, similar 
to the "universality classes" of critical phenomena [8] . 

One such distinction in two classes has been proposed by Angell [5], namely 
the distinction between "strong" and "fragile" glassformers. Plotting the logarithm 
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of the viscosity i](T) versus the normalized inverse temperature T g /T (the glass 
transition temperature T g is here defined somewhat arbitrarily from the condition 
r)(T = T g ) = 10 13 Poise), one observes that certain network-forming materials such 
as molten SiC>2 and molten GeC>2 simply follow straight lines, i.e. the temperature 
dependence of rj(T) can be described by an Arrhenius law, 

r)(T) = Tfco exp (Ejk B T) , (1) 

where r/^ is a constant and plays the role of an activation energy. Most other 
glassforming systems, however, in particular polymer melts, multicomponent metallic 
melts, and fluids formed from small organic molecules, behave differently. For these 
glassformers, the plot of log[^(T)] vs. T g /T is strongly curved. Following Angell [S], 
these systems are called "fragile glassformers" . 

There is ample evidence [3 [TU] that in fragile glassformers the initial stages of 
slowing down, when the structural relaxation times grow from the picosecond scale 
by several orders of magnitude, can be described rather well by mode coupling theory 
(MCT) [3j, although some aspects of this theory are still under discussion [TT] . and 
there is no consensus on the behavior near T g [6j [1~3 [13 [14] . For the case of silica, 
computer simulation studies [151 1161 117j have shown that the relaxation dynamics at 
high temperatures can be well described by MCT, whereas at low temperatures an 
Arrhenius behavior is observed, as seen in experiments (note that the high temperature 
regime is almost not accessible by experiments). This indicates that, at least on a 
qualitative level, the "strong glassformers" SiC>2 exhibits a similar behavior for the 
temperature dependence of transport coefficients and structural relaxation as typical 
"fragile glassformers". Now, the question arises whether this is also true for the 
other prototype of a "strong glassformer" , namely GeC>2. While molten silica has 
been studied extensively, both by various experimental techniques and by computer 

simulations [13 [13 [12 [13 [13 H [233 [23 [^ 

studies of molten and glassy germanium dioxide are less abundant, and this holds true 
for both experiments [351 IHHl 133 123 HOI HH H3 H3 HH HS1 HBJ and simulations 

03S3HSII53I5I1I53I53E1 1551- 
in the present work, we hence present a detailed Molecular Dynamics (MD) 
|56[ I57j study of molten and glassy Ge02 at zero pressure, using a pair potential 
model that has been recently proposed by Oeffner and Elliott [47 . In order to check 
whether the Oeffner-Elliott (OE) potential provides a chemically realistic modelling 
of GeC>2, we perform also "ab initio" Car-Parrinello Molecular Dynamics (CPMD) 
simulations 58, 59, 60] and compare various structural and dynamic quantities as 
obtained from classical MD using the OE potential with those from the CPMD 
calculations. Moreover, our simulation results are also validated by comparison to 
experimental data. 

In Sec. 2 we summarize the models and methods of the simulation, while Sec. 3 
is devoted to a description of the static properties of molten and glassy Ge02 (partial 
pair correlations and structure factors, ring statistics and angular distributions, 
etc.). Section 4 presents selected information on dynamic properties (mean square 
displacements, intermediate incoherent scattering functions), while Sec. 5 summarizes 
some conclusions. 



Molecular Dynamics study of molten and glassy germanium dioxide 



3 



2. Models and simulation methods 

2.1. Classical MD 

In a classical MD simulation, all degrees of freedom due to the electrons are 
disregarded, as well as quantum effects due to the ions (which need to be included 
for a correct description of thermal properties of glasses at temperatures far below 
the glass transition temperature). One simply solves Newton's equations of motion, 
which is conveniently done applying the velocity form of the Verlet algorithm J55J ISU • 
Forces are computed using the OE potential [4"7] , 

2 

VcfiiVij) = ^^+A af3 eM-B a /3r t3 )+C a f 3 r- j 6 a,(3e Ge,0 .(2) 

Here, = \fi — fj\ is the distance between a pair of particles at positions r, and 
fj. The first term on the right hand side of Eq. ([2J describes Coulomb interactions, 
with e the elementary charge and the values qq c = 1.5 and go — —0.75 for the 
partial charges of germanium and oxygen ions, respectively [47j . The second and 
third term in {2J form a Buckingham potential and describe the short-range part of 
the potential. The constants B a p and C Q( g are [47] Aq o = 208011. 52eV, -BgcO = 
6.129329A" 1 , C Go0 = 236.653eVA 6 , A QO = 7693.522eV, B 00 = 3.285108A- 1 , 
and Coo = 131. 09 eV A 6 . The Buckingham terms for the Ge-Ge interaction are 
set to zero. The OE potential was derived from quantum-chemical calculations of 
Ge04 tetrahedra, using also experimental data from the a-Ge02 crystal structure at 
T = 300 K as input information. We note that analogous procedures for the chemically 
similar case of Si02 have led to the potential due to van Beest, Kramer and van 
Santen ("BKS potential") .61!], which has proven useful to reproduce a great variety 
of experimental results rather accurately [T51 [22J 1251 [25J [23] . Note that for Si02 the 
effective charges are different (qsi — 2.4, qo — —1-2), despite the chemical similarity. 
Thus, a combination of the OE and the BKS potential would not be suitable for the 
description of oxide melts containing both Ge02 and Si02, since the 0-0 interaction 
is modelled differently in both cases. 

Several other potential models were proposed in the literature [511 E21 [63] , but 
structural properties of liquid Ge02 derived from these potentials are not in good 
agreement with experiment, and hence these potentials were not used in the present 
study. 

While the short-range part of Eq. (fT]) was cut off and shifted to zero at a distance 
r c = 7.5 A [5J], the Ion g-range Coulomb interactions were treated by Ewald summation 
methods [55J [53]. The equations of motion were solved for systems of N = 1152 
atoms using an integration time step of 1.23 fs and periodic boundary conditions in 
all three spatial directions. The simulations in the NpT ensemble at constant zero 
pressure, p = 0, yielded linear dimensions L(T) of the cubic simulation box in the 
range 26.6 A < L(T) < 28.4 A for temperatures in the range 2530 K < T < 6100 K. 
Pressure was kept constant using the Andersen barostat [66j . Constant temperature 
was realized by coupling the system periodically to a stochastic heat bath [55] . Note 
that the runs in the NpT ensemble were only used to create well equilibrated initial 
configurations for runs in the microcanonical NVE ensemble (V denoting the system 
volume and E its internal energy). Using force parallelization with message passing 
interface (MPI) routines, an efficient use of the Julich multiprocessor system (JUMP) 
with 32 processors used in parallel was possible. Equilibration times t e spanned the 
range from 48.9ps (40000 time steps) at T = 6100K to 11.97ns (almost 10 7 time 
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steps) at T = 2530 K, to generate 8 initial configurations, which then were propagated 
in the NVE ensemble for the same time interval t e , during which structural and 
dynamical properties were recorded. Note that the time t e was chosen such that the 
slower species (Ge) moved on average a distance of 5.5 A at each temperature. Further 
implementation details are documented in Ref. [64] . 

2.2. CPMD 

A important issue for the CPMD simulations of Ge02 is the quality of the pseudo 
potentials, which are a necessary input for the CPMD method [551 HH]- While 
we found that a pseudo potential due to Goedecker et al [67] was computationally 
too demanding for our purposes, a pseudo potential based on the general gradient 
approximation (GGA) with the BLYP exchange-correlation functional in the Troullier- 
Martins parametrization [68j was found to be satisfactory. As an energy cutoff for the 
plane waves E cni — 75 Ry was used, similar as in related work for Si02 t 28j. The time 
step was 0.0726 fs. For the thermostatting of the system, we used Nose- Hoover chains 
[69] for each ionic degree of freedom as well as for the electronic degrees of freedom to 
counterbalance the energy flow from ions to electrons [70] . The parameters used for 
the Nose-Hoover chains can be found in a previous publication [28 . 

An important problem in CPMD simulations of amorphous systems is the 
generation of suitable initial configurations. While in the case of Si02, it was found 
useful to start from classical MD simulations using the BKS potential [61] and relax 
these configurations to new equilibrium states by CPMD [28j [71] [72] , in the case of 
Ge02 (using the OE potential [47]) such a procedure did not converge [64]. The reason 
for this failure is that the differences between equilibrated atomic configurations using 
either classical MD or CPMD methods for Ge02 are slightly larger than for Si02, 
as far as interatomic distances, angles etc. are concerned. At the temperatures of 
interest (T = 3760 K and T = 3000 K), which are far above the melting temperature 
T m of Ge02 (T m = 1389 K [73]) it is also too time-consuming to start from a crystalline 
configuration and melt it in a CPMD run; thus we decided to start from configurations 
generated by classical MD at T = 7000 K, where subsequent equilibration by CPMD 
turned out to be feasible (for 60 particles this took 53000 CPMD steps, while for 
120 particles 21000 CPMD steps were sufficient, using periodic boundary conditions 
throughout). Then the temperature was lowered in a single step to T — 3760 K (for 
N = 60) or T = 3000 K (for N = 60 and N = 120), respectively. At T = 3760 K, runs 
over 171000 time steps for equilibration and production were performed corresponding 
to a real time of 12.4 ps. At T = 3000 K, we did runs over 340000 time steps for the 
system with 60 particles and 420000 time steps for the system with 120 particles, thus 
covering a time range of 24.7ps and 30.5 ps, respectively. In order to obtain better 
statistics, we averaged over 6 independent simulation runs for each system size and 
temperature considered. 

The density was chosen to be p — 3.45 g/cm 3 , similar to the equilibrium density 
resulting from the classical MD simulations in this temperature range, in order to be 
able to compare MD and CPMD results at essentially the same density. This choice 
implies linear dimensions of the simulation box of L — 10.023 A for N = 60 and 
L = 12.629 A for N = 120. Since the periodic boundary condition does significantly 
affect the structure and correlation functions for distances that exceed L/2, the 
smallness of N and L clearly is a major disadvantage of our implementation of CPMD, 
and prevents us from a meaningful study of intermediate range order by CPMD. 
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Figure 1. Density of Ge02 plotted vs. temperature. Well-equilibrated MD 
results (diamonds), using the OE potential, are shown in the temperature range 
6100 K> T >2530K. The dotted line connecting the data points serves only as a 
guide to the eye. The MD data shown for T < 2750 K result from cooling runs with 
two different cooling rates, using well-equilibrated configurations at T = 2530 K 
as a starting point [cf. Eq. ((3}] . All the simulation results were obtained at zero 
pressure. Experimental data from Riebling | 35 | and Dingwell et al I38| are shown 
for comparison. 



Application of novel versions of ab initio MD, suitable to simulate significantly larger 
systems [73], is desirable, but must be left to future work. 

Finally, we mention that sometimes the generated configurations had to be 
discarded "by hand", when they contained well-identifiable O2 molecules disjunct 
from the remaining germanium oxide network (which then necessarily has coordination 
defects, of course). It is clear that at T — 7000 K such chemical disintegration of GeC>2 
may be a physically meaningful effect. But we are interested in the properties of GeC>2 
at lower temperatures, where these separate O2 molecules should no longer occur, but 
rather should be integrated into the network structure again. Even at temperatures 
around 3760 K, these O2 molecules are expected to be unphysical. If at much lower 
temperature we were able to equilibrate the samples over many orders of magnitude in 
time, we should see automatically that such defects anneal out again. However, since 
this is practically impossible, the somewhat biased selection of physically meaningful 
configurations had to be made. 

With respect to other implementation details, we closely followed the procedures 
of Bcnoit et al [28] (see also [64]). We only note that, in our case, the CPU time 
required for the CPMD is a factor of 358000 higher than that needed for the classical 
MD, using the same multiprocessor system and the same system size for both methods 
[51] , Therefore, only a rather restrictive use of CPMD was feasible. 
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3. Static properties of molten and glassy Ge0 2 

As discussed in Sec. 2.1, equilibration was done in the framework of classical MD 
using the NpT ensemble which allows to record the temperature dependence of the 
density (Fig. [I}. In our MD simulation, the lowest temperature which could still be 
equilibrated with manageable effort was T — 2530 K. This temperature corresponds to 
almost twice the melting temperature [73] , while experimental data are only available 
at much lower temperatures. Therefore, we used states at T = 2750 K for further 
cooling down the samples (note that the states at T = 2530 K were not available 
yet when these cooling runs were performed). To this end, temperature was linearly 
decreased according to 

T{t) = 2750 K- Qt, (3) 

with cooling rates Q = 2.25 x 10 13 K/s and 1.13 x 10 12 K/s. As in the case of Si0 2 
[T5I |2"2"] , the cooling rates available in MD exceed those of the experiment by many 
orders of magnitude, and a meaningful extrapolation to these very small experimental 
cooling rates is not possible. Although the presence of a density maximum (as is known 
to occur in SiC>2 [75]) somewhere around T — 2000 K cannot be excluded, it seems 
very unlikely that for slow cooling rates the simulated densities for T < 1700 K would 
decrease enough to match the experimental data. So we attribute the larger part of the 
mismatch between simulated and experimental melt densities to the inadequacy of the 
OE potential to predict the density very accurately! However, such a 5% discrepancy 
in the density is not uncommon when classical pair potentials are used. 

Surprisingly, at T = 300 K the experimental density is p cxp « 3.65 g/cm 3 and 
the simulated one (with the slowest of our cooling rates) p s j m w 3.70 g/cm 3 , thus only 
1.37% higher. However, this good agreement presumably is due to a lucky cancellation 
of errors (freezing in a too high density due to the inaccurate potential, partially 
compensates for not reproducing the rapid variation of the density of supercooled 
GeC>2 around T — 1000 K due to our by far too fast cooling) . This example again shows 
that fits or misfits of isolated experimental data points by simulations are unsuitable 
to judge the quality of potentials and/or simulation procedures. 

A more detailed information on the static structure, also available via neutron 
scattering experiments, is the static structure factor. Since we deal here with two 
species, it is appropriate to consider partial structure factors S a p(q) (a, (3 = Ge, O) 

N a N/3 

S a p{q) = jz ( £ exp(iq- ry)} • (4) 
i=i j=i 

Note that fluids and glasses are isotropic and hence S a p(q) depends only on the 
absolute value q = \q\ and not on the direction of the scattering vector q. Using 
suitable isotopes, all partial structure factors for GeC>2 have recently been measured by 
Salmon et al [45 , 46 . Figure [2] reveals a very good agreement between our simulation 
results and these data. 

Standard neutron scattering yields a scattering intensity weighted with the 
scattering lengths b a , bp as follows 

S(q) = r „, 2 2^ Kbp S al3 {q) . (5) 

Oi Qi/3e{GOi0} 

Using 76J &Gc = 8.185 fm, bo — 5.803 fm one can compute from Eqs. (HKHJ) the neutron 
scattering structure factor from the simulation and compare it to corresponding 
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Figure 2. Partial neutron scattering structure factors S a p(q) plotted 
vs. wavenumber q, comparing the present MD simulation to the experimental 
data of Salmon et al [25] Hg] at T = 300 K. 






Figure 3. Partial structure factors <SGeGc(<?)> part a), Sq c o(q)^ part b), and 
>Soo(<?)i part c), plotted vs. q and four temperatures, as indicated. 
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experimental data [32] without any adjustable parameters whatsoever. Also for this 
comparison [5U1 ES] the general agreement between simulation and experiment is 
rather good; both predict a "first sharp diffraction peak" (FSDP) [JJ H3 [H] a * about 
<Zmax ~ 1-55 A~\ which can be attributed in real space to the linear dimension of two 
GeC>4 tetrahedra sharing a corner (see below), £ = 27r/<7 max w 4.05 A. 

When we compare to SiC>2 [13 [78] we note that in SiC>2 the FSDP occurs at a 
slightly larger value, g m ax ~ 1.7 A -1 , implying a somewhat smaller linear dimension 
of the two corner-sharing SiC>4 tetrahedra (note that the "chemical rules" [JJ for 
the formation of perfect binary continuous random networks, with a cation in the 
center of a tetrahedron and oxygens at the corners, such that each oxygen is shared 
by two neighboring tetrahedra, are identical for SiC>2 and GeC>2, of course). But 
a more interesting difference is the fact that SiC>2 shows a second well-developed 
peak, at about q' max ~ 3 A -1 , which corresponds to a peak in GeC>2 at about 
9max ~ 2.6 A -1 . While in the total neutron scattering structure factor this peak 
is hardly distinguishable from the noise, the partial static structure factors (Figs. [21 
[3]) reveal that actually this is the main peak in the structure, corresponding to a 
distance £' = 2n/q' may , w 2.4 A. This distance, however, cannot be attributed to any 
interatomic distance in the structure of GeC>2 directly. It rather corresponds to the 
period of the oscillatory decay of the partial pair correlation functions g a p(r) m rea l 
space at large distances (Fig. 0]). These correlations are obtained from the simulated 
configurations from their definition 

M0 = ^EEl3 J ( r -| f '- f 4 a,/3-{Ge,0}, (6) 
i=l j=i 7rr 

where N a ,p = V/{N a N ) if a ^ while M aa = V/[Af(N a -i)], V being the volume of 
the simulation box. The correlation functions g a 0{r) and S a p(q) are related via [6l [79] 

S aP (q) = 1 + (N/V) J [g aP (r) - 1] exp{iq- r) dr. (7) 

In the following, we shall focus on g a p(r) rather than on S a p(q). Ref. [44] did give 
some estimates of the nearest neighbor distances of the various types of pairs, which 
are included in Fig. [U indicating a reasonable agreement with the simulation. Note 
that the Ge-0 distance (about 1.73 A) clearly is the smallest distance occurring in 
the structure, and the sharpness of this peak [note the ordinate scale of Fig. |4}d) in 
comparison to that of Fig. [4^,!] reveals that the Ge04 tetrahedra are fairly rigid. Only 
for the Ge-Ge distance a slight systematic discrepancy between MD and experiment 
is visible. Comparing to the CPMD results (Fig.O, however, this discrepancy seems 
to be removed. 

It is also possible to compare CPMD results for S a p(q) and S(q) with the 
corresponding MD results [64] : these comparisons strengthen the conclusion that one 
can also draw from Fig. [5[ namely that MD yields a rather accurate description 
of the local structure of molten Ge02- Note that slight discrepancies in goo(r) 
for r > 5 A should not be taken very seriously, because at these distances CPMD 
suffers from finite size effects, as noted above. More interesting is the difference 
(emphasized in the insert of Fig. concerning the feature near 1.5 A. Testing carefully 
different equilibration times it was possible to show that too short equilibration of 
CPMD yields such prepeak in goo( r ) which is too high rather than too low [64] . 
Therefore, this difference between the CPMD and the MD results is probably a 
real effect, at least it is not an artifact of too short equilibration. Of course, one 
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Figure 4. Partial pair correlation functions g a p(r) for Ge02 plotted vs. r for 
various temperatures, as obtained from the classical MD results. The broken 
vertical straight line indicates the estimates of nearest neighbor distances of the 
Ge-Ge- pairs (a), Ge-O-pairs (b) and O-O pairs (c) as extracted experimentally 
from measurements of partial structure factors at T = 300 K |44| . 
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Figure 6. Schematic picture of a ring of length n = 6, illustrating also the 
definition of an angle 8q c q e q c and the distance d between neighboring Ge atoms 
(a). A ring of length n = 2 and the angle £*GeOGc is sketched in (b), and the 
tetrahedral angle #oGeO m ( c )- 



can question the accuracy of CPMD somewhat on other grounds: other "ab initio" 
studies of the GeC>2 structure (50] [80] using different pseudopotentials and system 
preparation procedures predicted somewhat different results (e.g. the Ge-0 distance 
r GcO — 1-69 A [80 or rccO = 1-78 A [5U], while we obtain 1.71 A and the experimental 
value is 1.73 ± 0.03 A [SS gj). 

We now turn our attention to the analysis of structural features on intermediate 
length scales. To this end, we recall the concept of "ring statistics" [6] [22]. One 
considers the shortest closed paths in the network of covalent bonds, starting from an 
oxygen atom (Fig. [5]). The length n of a ring is then the number of cations (Ge in the 
present case, or Si in the case of silica [22]) that one passes before one returns to the 
starting point. Figure [6] shows, as an example, n = 6 (left) and n — 2 (middle part). In 
Si02, it has been found that the angles between atoms in a ring with n = 2 and n = 3 
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Figure 7. Probability P(n) that a ring of length n occurs, plotted vs. n, for 
T = 3000 K and T = 3760 K, comparing MD and CPMD (a), and the same 
comparison including an MD study where N = 60, as in the CPMD calculation 
(b). 



differ appreciably between the classical MD simulation and its CPMD counterpart 
[ZU E2] , and this affects also significantly the probability P(n) that a ring of length n 
occurs in the structure (in thermal equilibrium) . At first sight one might conclude that 
a similar effect occurs for GeC>2 , too (Fig. [7^) , but a closer analysis reveals that most 
of the differences between CPMD and MD stem from the fact that the former suffers 
from finite size effects (Fig. Wp)'- when we use N = 60 in the MD calculation, we find 
almost perfect agreement with the CPMD calculation that uses N — 60 as well. Also 
the strong difference between the CPMD results for ^=60 and iV=120 show that one 
cannot trust the CPMD results for P(n), due to these dominating finite size effects. 
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Figure 8. Distribution functions P(8) of various angles 8, obtained from MD for 
a wide range of temperatures, as indicated. Case a) shows the Ge-Ge-Ge angle, 
case b) the Ge-O angle, case c) the O-Ge angle and case d) the O-O-O angle. 



Clearly, for a quantity that depends sensitively on the order on intermediate length 
scales like P(n) it is more important to choose a large enough system rather than to 
work with very realistic descriptions of the forces, as provided by CPMD. 

While CPMD hence is less useful for the study of properties that depend 
sensitively on medium range order, it clearly is of great interest for the assessment 
of local properties, such as the distributions of angles between the "bonds" in the 
structure. These distributions have also been obtained by MD for a wide range of 
temperatures (Fig. [8]). The definition of the Ge-Ge-Ge angle is indicated in Fig. [6^,); 
other angles are defined analogously. A remarkable feature is that all distributions, 
with the exception of the tetrahedral angle O-Ge-O, have a double peak shape, and are 
rather broad. Only the distribution of the tetrahedral angle tends towards a Gaussian 
shape, as the temperature is lowered, and gets somewhat sharper; in a random network 
structure formed by ideal tetrahedra only, this distribution would be a delta function, 

6(0 - 0tetr) With tetr = 109°. 

The relative weight of the peak at 8 = 60° of the Ge-Ge-Ge angle decreases with 
decreasing temperature, as well as the weight of the peak at 8 = 90° for the Ge-O-Ge 
angle distribution. A consideration of the geometry of the rings (Fig. [6]) immediately 
shows that the peak of P(8) for the Ge-Ge-Ge angle can be attributed to rings with 
n = 3, and similarly the peak of P(8) for the Ge-O-Ge angle at 8 = 90° is due to 
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rings with n = 2. Such small rings can be frequently observed in the structure of 
Ge02 at high temperatures, while at low temperatures the network becomes much 
more regular, and the density of all small rings decreases significantly. 

For T — 2530 K, the position of the main peak of the distribution P(9) for the 
Ge-O-Ge angle is 133°. It is gratifying that this number coincides with corresponding 
experimental estimates (39) [81] ■ This agreement is a further indication that the OE 
potential is able to provide a rather realistic description of the structure. 

The side peaks of Fig. [5K,b, tend to disappear at the physically relevant 
temperatures, i.e. the number of rings with n — 2 and n = 3 becomes significantly 
smaller with decreasing temperature. This is also indicated by the temperature 
dependence of the ring length distribution P(n) [64]. The main peaks of the GeGeGe 
and GeOGe distributions seem to stay rather broad, as expected due to the disorder 
in the network structure. Only in the various crystal structures of Ge02 at low 
temperatures we would expect very sharp distributions of all angles; in the glass 
structure only the distributions of the angles inside a tetrahedron become rather sharp 
at low temperatures. 

In this respect, the distribution of the angle 9 between 0-0-0 bonds is special: 
Fig. EJi shows that there two peaks occur, which clearly persist at low temperatures. 
The obvious explanation is that there are two distinct possibilities: the peak at 9 = 60° 
can be attributed to oxygen atoms belonging to the same tetrahedron, while the peak 
at 9 <w 110° is due to oxygens belonging to two neighboring tetrahedra. In fact, as 
temperature decreases the structure of a single tetrahedron approaches more and more 
that of an ideal tetrahedron, whose faces are perfect triangles, having angles of 60°. 
In view of this, the observation that the peak at 9 = 60° becomes clearly sharper with 
decreasing temperature is not surprising. 

Now we turn to the comparison of these angular distributions to the corresponding 
CPMD predictions (Fig. [9]). The general shape of these distributions is very similar, 
with the exception of the Ge-O-Ge angle, where the side peak at 90° (due to rings 
with n = 2) is broadened into a shoulder only, indicating that the OE potential 
overestimates in particular the rigidity of this structural element (a schematic picture 
of a ring with n — 2 is shown in Fig.lBJ})). Wc also note that the CPMD distributions 
are always somewhat broader than the MD results at the corresponding temperature. 
This indicates that the CPMD calculation, if we could parametrize it in terms of an 
effective pair potentials having the OE or BKS form, would yield a systematically 
softer potential. In fact, if one compares the CPMD calculation at T = 3000 K to 
the classical calculation at T = 3760 K, the differences are much smaller [53 . Of 
course, we do not wish to imply that the differences between CPMD and MD could be 
fully eliminated by a renormalization of the temperature scale: for the main peak of 
the Ge-O-Ge angle distribution, CPMD at T = 3000 K implies a peak at about 129°, 
while the MD calculation yields a peak at about 133° (this value depends much less on 
temperature than the CPMD peak position does). We have also done MD simulations 
with N — 60 particles only, to rule out that the differences seen in Fig. [9] simply are 
due to finite size effects [64]. Figure [9^, also indicates for the Ge-Ge-Ge distribution 
that CPMD result for N = 60 is only slightly different from that at N = 120 (for the 
other distributions the differences arc even smaller). 

As a conclusion of this section we may state that the OE potential predicts slightly 
too rigid structures in comparison to CPMD, and this difference is most pronounced 
at rather high temperatures. However, the overall agreement between the structure as 
predicted by the OE potential and the structure resulting from CPMD is very good. 





Figure 9. Comparison between MD and CPMD results at T = 3000 K for the 
distribution functions of various angles, Ge-Ge-Ge (a), Ge-O-Ge (b), O-Ge-O (c) 
and O-O-O (d). All MD results refer to TV = 1152, while the CPMD are results 
are for N = 120 (for the Ge-Ge-Ge distribution also the CPMD result for N = 60 
is shown). 



The same conclusion emerges also from an analysis of the distribution of coordination 
numbers [64] • The comparison to experimental data, whenever available, also suggests 
the statement that the OE potential provides a reasonably accurate description of the 
static structure of molten and glassy GeC>2. 



4. Dynamic properties of GeC>2 melts 



From the MD runs in the NVE ensemble, it is straightforward to record both the mean 
square displacements (MSD) of a tagged particle of type a {a — {Ge, O}) [41161179]. 



rl(t) 



£(!*(*) -*(o)l 



i=l 



and the intermediate incoherent scattering function 



F s a (l,t) = ±- £ (exp-Hg- [n{t) - f-(0)]} 



(8) 



(9) 



The MSD allows to estimate the self-diffusion constants, applying the Einstein relation 



D a = lim 

t — >oo 



(10) 
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Figure 10. Log-log plot of the MSD for Ge(a) and O(b) versus time, for 
temperatures ranging from T = 2530 K to T = 6100 K. 
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Figure [10] shows our MD data for the MSD. One sees the standard behavior, 
familiar from MD simulations for SiC>2 [15j and many other systems [6]. At very 
short times, a ballistic regime is seen ((r^(t)) oc i 2 ). Then, at high temperatures, 
a rapid crossover to the linear diffusive regime occurs ((r 2 (t)) = 6D a t), while at 
lower temperatures, i.e. in the range 2530 K< T < 3250 K, a plateau is observed at 
intermediate times, where the MSD does not increase, but rather stays constant at 
about (r 2 (t)) « 0.5 A 2 . This plateau commonly is interpreted as the onset of the "cage 
effect" [5J |S] : each atom sits in a "cage" formed by its nearest neighbors and the lower 
the temperature the more time it takes until the atom can "escape from the cage" . 
Of course, such mobility implies that the network of bonds in the random network 
structure is not rigid, sometimes a bond "breaks" 15J and coordination defects appear, 
which later can anneal again. 

The MSD, as obtained from MD simulation with the OE model, can be also 
compared to corresponding CPMD results. Figure [TT] shows a behavior which is not 
surprising at all, in view of our findings for static properties as described in detail 
in the previous section: The time dependence of the MSD found for T = 3000 K by 
CPMD superimposes almost exactly with the MD results for T = 3760 K, reflecting 
again the finding that CPMD is essentially equivalent to the use of pair potentials 
that are slightly softer than the OE potential but otherwise very similar. As a further 
caveat we mention the effect of the Nose- Hoover thermostat (needed in CPMD, not 
in MD), which may have speeded up slightly the CPMD dynamics, though we do not 
have any real evidence that this effect is already important on time scales up to 20 ps 
that are shown in Fig. [TT] 
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Figure 12. Self-diffusion constants of Ge and O in Ge02 melts plotted vs. inverse 
temperature. For comparison, results for Si and O in Si02 melts (taken from 1151 ) 
are included. Straight lines indicate fits to the Arrhenius relation, Eq. mil l. Also 
an Arrhenius fit resulting from experimental viscosity data 1351 via Eq. 1121 is 
included. 



In Fig. [T^] a plot of the diffusion constants is presented, choosing a logarithmic 
ordinate scale and inverse temperature as abscissa, so Arrhenius relations show up via 
straight lines, since then [compare to Eq. ([T])] 

D a = D a ^ exp[-£ a , Q /(fc s T)] . (11) 

The activation energies resulting from the fits in Fig. fTS] are -E a .Gc = 3.41 eV and 
E a .o = 3.25 eV. As can be also infered from Fig. [T31 oxygen diffuses slightly faster 
than Ge, and this difference becomes slightly more pronounced with decreasing 
temperature, due to the slightly higher activation energy of Ge. A similar behavior is 
well-known for Si02 [TS] . 

While in the case of Si02 experimental data for self-diffusion constants £>si, 
Do are available, we are not aware of suitable data for Ge02- However, when we 
disregard the small difference between Do e and Do, a rough estimation of these 
diffusion constants is possible with the well-known Stokes- Einstein relation (6] [79] , 

D = kgT /\cKrjR\ , (12) 

with the constant c = 4 if one assumes slip boundary conditions for the particle 
diffusion in the fluid, r\ is the shear viscosity of the fluid, and R the radius of the 
diffusing particle. In principle, Eq. (|12[) is a result from hydrodynamics, and makes 
sense only if R is much larger than interatomic distances. However, in the spirit of the 
finding that often descriptions based on hydrodynamics work down to the molecular 
scale (see [55] for a recent example), Eq. (fT2"j) is used also for diffusing atoms or 
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Figure 13. Incoherent intermediate scattering functions F^ c (q,t) (a) and 
Fp(q,t) (b) for Ge02 at q = 1.55 A -1 , where the static structure factor S(q) 
exhibits the first sharp diffraction peak, plotted vs. time (on logarithmic scale) 
for a broad range of temperatures (the leftmost curve corresponds to T = 6100 K, 
the rightmost curve to T = 2530 K, temperatures in between are 5200 K, 4700 K, 
4300 K, 4000 K, 3760 K, 3580 K, 3400 K, 3250 K, 3100 K, 3000 K, 2900 K, 2750 K, 
and 2640 K). 
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molecules. Using then for R the Ge-0 nearest neighbor distance, R = 1.75 A, the 
experimental viscosity data of Riebling [35] are readily converted into the self-diffusion 
constant, and the resulting Arrhenius fit (implying E a = 3.565 eV) is also included in 
Fig. 1121 and in very good agreement with our simulations. 

In contrast to our results, previous simulations [5T][52] gave activation energies in 
the range between 1 eV and 1.2 eV. There are many indications that the potential used 
by Hoang [51] cannot describe GeC>2 as accurately as the OE potential does; moreover 
his results presumably suffer from aging effects due to insufficient equilibration. The 
latter criticism also applies to the study of Micoulaut et al [52] . where the system 
configurations were taken from one cooling run applying a cooling rate of 2.5 x 10 12 K /s, 
although states in the temperature range from 940 K < T < 2480 K were considered. 
It is clear that such configurations are far from equilibrium, even at T = 2940 K 
these data (52] do not show any sign of the cage effect, and at the lower temperature 
(T = 940 K), which is only about 100 K higher than the experimental glass transition 
temperature, the structural relaxation time is only of the order of ns, which proves 
that the melt is in a state very far from equilibrium, in the initial stages of aging. 

We now turn to the analysis of intermediate scattering functions (Fig. [T3|) . Again 
we note the qualitative similarity of these curves to data for many other glassforming 
fluids [Hd!]. While at high temperatures the decay of F"(q,t) resembles a simple 
exponential, for T < 3400 K the decay occurs in two steps, due to the cage effect. The 
so-called "/3-relaxation" is the time regime around the plateau, while the final decay 
from this plateau to zero is called "a- relaxation" [31 [6]. Whereas, at high temperatures 
F"(q,t) decays to zero on the ps timescale, at the lowest accessible temperatures the 
"lifetime" of the plateau extends into the ns time range already. In order to define the 
structural relaxation time T a (q, t), we follow Ref. [T7] by requesting that for t = T a (q, t) 
the scattering function has decayed to a value of 0.1. Thus, 

F s a (q,t = r a (T)) = 0.1, a = Ge,0; (13) 

here we have omitted the argument of the structural relaxation time, since in the 
present context only the value of q, corresponding to the location of the FSDP in the 
static structure factor, is of interest. 

Figure [14] presents a log- log plot of tq c (T) and to(T) against T — T c , T c = 
2490K±100K being an estimate for the mode coupling [3] critical temperature. Of 
course, it is well-known that no real divergence of r a (T), implying an ergodic-to- 
nonergodic transition and a divergence of the viscosity at T c , can occur in real 
glassforming fluids 6 : rather idealized mode coupling theory [3] is a kind of mean field 
theory for the dynamic correlation functions of glassforming fluids, which supposedly 
holds for T > T c but not too close to T c , since the predicted divergence at T c is 
rounded off. The standard albeit heuristic interpretation is that an infinite lifetime of 
the cage "imprisoning" of the particles is prevented by thermally activated processes, 
so-called hopping processes, which break up the cage even for T = T c and T < T c . So 
T c only plays the role of a crossover temperature, where the temperature dependence 
of T a (T) crosses over from the power law T a (T) cx (T — T c )~ 7 to an Arrhenius law. 
Note that we find the value 7 = 2.5 for the critical exponent. This value is slightly 
higher than the one estimated from simulations of silica [T7]. The insert of Fig. [TJ] 
shows that indeed at lower temperatures our estimates for r a (T) are consistent with 
an Arrhenius law. The activation energy is E a = 3.57 eV in this case, i.e. slightly 
higher than those determined for the self-diffusion constants (see Fig. IT2"]) . Of course, 
this crossover from a power law to thermally activated behaviour is by no means sharp, 



Molecular Dynamics study of molten and glassy germanium dioxide 



22 




1/7(10" 4 K" 1 ) 



10 1 10 2 10 3 

T-T C (K) 

Figure 14. Log-log plot of TQ e and to(T) in fluid Ge02 versus T — T c , using 
T c = 2490 K. The bold lines are fits with power laws oc (T - T C )T using 7 = 2.5. 
The insert shows the same data as log[r a (T)] versus l/T. The straight line in the 
insert is a fit with an Arrhenius law (oc exp[Ea/ (kgT)]) with E & = 3.57 eV. 



but actually rather gradual, and this smooth behaviour near T c necessarily prevents 
us from an accurate estimation of T c : data near T c may deviate from the straight 
line on the log-log plot due to the onset of the crossover, even if the estimate for T c 
is correct; thus, the precise range of temperatures for which T a (T) should be fitted 
to the power law is somewhat uncertain, and this leads to a considerable uncertainty 
about T c . 

As further evidence for the applicability of mode coupling theory to describe 
the dynamics of molten GeC>2 at high temperatures, Fig. Lt5l shows a test of the time 
temperature superposition principle [3, 6]. Of course, the /3-relaxation is not supposed 
to follow this a-scaling and thus the upper part of the curves splay out, while the decay 
of the plateau (for temperatures where a plateau exist) follow this scaling nicely. The 
quality with which this scaling holds clearly is comparable to that observed for silica 
and fragile glassformers. 

In SiC>2, Saika-Voivod et al |30j have tried to link the relaxation dynamics to 
structural anomalies, such as the occurrence of a density maximum [30j . These authors 
have found some evidence for the occurrence of a kind of liquid-liquid phase transition 
in fluid SiC>2 at suitable conditions of temperature and pressure, i.e. there should exist 
two phases of fluid SiC>2 with different densities and different structure of the random 
network of covalent bonds. If this interpretation is correct, Figs. [12] and [14] would 
suggest that one should seek a similar interpretation in molten Ge02, too. However, 
our data (and the experimental data) for the density of GeC>2 do not give any hint for 




Figure 15. Incoherent intermediate scattering functions " (</, t) plotted vs. the 
scaled time t/r a for Ge (a) and O (b). The temperatures shown are the same as 
those of Fig. [13] 
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Figure 16. Plot of the product T a (T)D a (T) vs. temperature, for a = Ge and 
a = O, as indicated. Simulation results (data points) are connected by lines to 
guide the eye. 



structural anomalies of molten GeC>2 similar to those of SiC>2 (see Fig. [I). 

Finally, we address the question to what extent the Stokes-Einstein relation is 
valid for germania. In Eq. (|12p . this relation was used to link experimental viscosity 
data [35] to the self-diffusion constants (Fig. [T2]). While it would be possible to 
estimate the shear viscosity from the time-correlation of the off-diagonal pressure 
tensor components [SSJ [57], and such an approach has been shown feasible also for 
molten SiC>2 [T5], the statistical errors of the resulting estimates are very large. Since 
one usually associates the shear viscosity with the structural relaxation time, we 
tentatively tested to what extent the product T a (T)D a (T) is constant (Fig. \W\i . 
Similar as in Si0 2 , where the product r](T)D a (T)/T could be studied [TSJ, one 
sees some increase of this Stokes-Einstein ratio, but the increase again is not really 
dramatic. Much more dramatic violations of the Stokes-Einstein relation have been 
observed for various materials near T g , and this has found a lot of attention in 
the literature (see e.g. [21 [531 1EH US]). All of our data are far above the melting 
temperature of GeC>2, and we can say that in this temperature regime a dramatic 
breakdown of the Stokes-Einstein relation does not occur. If a breakdown occurs 
near T g , this would imply that the good agreement between the experimental data in 
Fig- HH and our estimates for the self-diffusion constants is merely accidental. 

5. Conclusions 

In the present work, we have described the results of a simulation study of fluid GeC>2 
at zero pressure, based on extensive MD runs using the Oeffner-Elliott potential. In the 
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temperature region T > 2530 K the melt has been carefully equilibrated, while glassy 
structures of amorphous Ge0 2 at room temperature were also produced, cooling down 
the system with two different cooling rates which did produce only minor structural 
differences, however. To validate our potential, also an "ab initio" Car-Parrinello 
Molecular Dynamics (CPMD) study was performed, which is limited to even higher 
temperatures (T > 3000 K) and very small systems (N < 120 atoms). When we 
consider properties at small enough scales where systematic errors due to finite size 
do not matter, we find very good agreement between the MD and CPMD descriptions, 
for both static and dynamics properties. The most important distinction is that the 
effective potential to which CPMD corresponds is slightly softer than the OE potential. 
However, the OE potential is clearly more accurate than other (empirical) potentials 
that were used in the literature, and our results agree also rather well with the (albeit 
somewhat scarce) experimental data which are available so far. 

Having validated the OE potential, it would be interesting to use it both for a 
more complete study of the dynamics of molten and glassy Ge02, and for a careful 
study of Ge02 under pressure. This must be left to future work. Also it would be 
interesting to use CPMD to construct a new effective potential which is even more 
accurate than the OE potential. Such attempts have been made [64], but so far were 
not successful. 

Our results suggest that for T > 2500 K a crossover sets in for the dynamical 
properties from a thermally activated behaviour of various quantities to a behaviour 
described by mode coupling theory, similar to previous findings for silicon dioxide. In 
view of our results it is not too surprising that in the temperature range T < 1600 K no 
experimental evidence for mode coupling effects could be found [4T] . Thus, it would 
be very useful if measurements could be extended to higher temperatures. Also a 
measurement of self-diffusion constants would be useful, to test the finding that only 
a rather weak violation of the Stokes-Einstein relation occurs in Ge02 . 

Ge02 and Si02 are the two archetypical examples for strong glassformers. Our 
results imply that they behave qualitatively similar, including the crossover to mode 
coupling type behavior at temperatures far above melting. The latter crossover is 
also characteristic for fragile glassformers. However, it is interesting that, in contrast 
to typical fragile glassformers, the critical mode coupling temperature is above the 
melting temperature. It has to be seen in future studies whether the relative location of 
critical temperature and melting temperature can be used as a criterion to distinguish 
fragile and strong glassformers. 

In conclusion, we hope that the present work helps to sort out the many questions 
concerning the possible universality (or lack thereof) in glassforming fluids, and 
stimulate further experimental and theoretical work on the above issues. 
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